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We study the dynamics of simple reactions where the chemical species are confined 
on a general, time-modulated surface, and subjected to externally-imposed stirring. 
The study of these inhomogeneous effects requires a model based on a reaction- 
advection-diffusion equation, which we derive. We use homogenization methods to 
show that up to second order in a small scaling parameter, the modulation effects 
on the concentration field are asymptotically equivalent for systems with or without 
stirring. This justifies our consideration of the simpler reaction-diffusion model, 
where we find that by modulating the substrate, we can modify the reaction rate, the 
total yield from the reaction, and the speed of front propagation. These observations 
are confirmed in three numerical case studies involving the autocatalytic and bistable 
reactions on the torus and a sinusoidally-modulated substrate. 
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1. INTRODUCTION 

We investigate the dynamics of the logistic and bistable reactions on a generalized, time- 
varying substrate with stirring. Such spatially-inhomogeneous problems call for the solution 
of a reaction-advection-diffusion equation, and much chemical and biological activity in fluid 
flow can be modelled by such equations. In particular, problems concerning autocatalytic 
chemical reactions [I] and population dynamics [2", "3] possess a logistic growth function as 
a reaction term, and thus satisfy a Fisher-KPP type of equation. Other, more complicated 
growth functions can be used to model a variety of phenomena, including the spread of 
insect populations, or the propagation of electro-chemical waves in organisms [SJ. The 
original motivation of this study was to understand the effects of wave modulation on the 
population dynamics of plankton, although the language and techniques we use are quite 
general. 

Before deriving and analyzing our model, we place our work in context by examining 
several streams of work that are relevant. It is known that a growing domain can modify 
biological pattern formation, as evidenced by the work of Newman and Frisch |4j. This has 
given impetus to the study of react ion- diffusion equations on growing domains [3 |6] in one 
dimension. Logically, this has led to the study of such problems on manifolds embedded in 
three dimensions. In multiple dimensions, one considers either the effect of curvature |7l El [9] , 
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or the twin effects of domain growth and curvature [TOt [TT] . The paper of Plaza et al. [10] 
is particularly relevant to the present work. In it, the authors derive the reaction-diffusion 
equation for a class of manifolds, and then study pattern formation on growing domains. 
We re-work their derivation to include the most general two-dimensional (differentiable) 
manifold possible, and then shift the focus from pattern formation to reactions in the pres- 
ence of stirring. The geometric formalism of Aris [I2] is central to our derivation. These 
references [101 [H] consider the 'geometric sink', that is, the notion that a growing domain 
can act as a sink for the chemical reaction. We extend this idea to growing and modulating 
domain sizes and examine the effects of the sink through numerical simulation. 

The notion of flow-driven reactions is not new. In [13], the effects of chaotic advection 
on the Fitz Hugh-Nagumo model are considered. The flow is found to produce a coherent 
global excitation of the system, for a certain range of stirring rates. The effects of flow can 
also induce distinctive spatial structure in the chemical concentration; this is studied in [13] . 
The same authors examine the single-component logistic or Fisher-KPP model in [13] . There 
the focus is on regime-change, namely how the rate of chaotic advection affects the spatial 
structure of the concentration. For slow stirring / fast reactions, a spatially inhomogeneous 
perturbation decays rapidly, and the equilibrium state is reached rapidly. On the other 
hand, for fast stirring / slow reactions, the perturbation persists, and a filament structure 
propagates throughout the domain. Nevertheless, the asymptotic state is still a stable 
homogeneous one. This system is simplified and the the transition treated as a bifurcation 
problem in [T3]. The focus of these papers is on local temporal and spatial structure. It is 
therefore salutary to examine the paper of Birch et al. [TB] , where the authors examine the 
averaged effect of a non-constant growth rate on the dynamics of the stirred Fisher-KPP 
equation. The authors use the theory of estimates to obtain bounds on the reaction yield, as 
a function of the stirring and the non-constant growth rate. When the mean growth rate is 
negative, the previously-unstable zero state of the Fisher-KPP equation can become stable; 
then the catalyst fails to propagate the reaction. In the Birch paper, the inhomogeneous 
growth rate is the consequence of an inhomogeneous distribution of nutrient in a plankton 
population. It could, however, be the result of placing the population or chemical species on 
a modulating surface, which is the subject of our report. Indeed, our results demonstrate 
the possibility of increasing the reaction yield by surface modulation. 

If the flow field or the modulation have small length scales compared to the domain 
size, then homogenization theory naturally presents itself as a tool for understanding the 
effects of flow and modulation in an averaged sense pJTj. The small scales are bundled up 
into an effective diffusion constant, and the model reduces to a more manageable equation 
involving a diffusion operator. Such an approach has been used in the linear advection- 
diffusion equation [TH [191 120], for linear reaction-advection-diffusion equations [21], and 
for non-linear equations We propose and justify the application of this theory to a 

linearized, advective Fisher-KPP equation on a time-varying manifold, and compute the 
effective diffusivity for the torus. The effective diffusivity for a manifold is, in general, a 
function of position, although in the fortuitous case we consider, this dependence vanishes. 
As in [23], the asymptotic state of the system is homogeneous, and the purpose of our 
homogenization calculation is therefore to study the speed at which this state is reached. 

This paper is organized as follows. In Sec. |2] we discuss the autocatalytic reaction in 
the homogeneous case, and then generalize this formalism to consider reactions on time- 
dependent, spatially-inhomogeneous surfaces. In Sec. [3] we derive conditions on the metric 
tensor for homogeneous solutions to exist on a general surface. We then outline a separation- 
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of-scales technique that enables us to compute the spatial distribution of concentration as the 
solution of a diffusion equation. We discuss the case of shear flow on a modulating torus, 
where the equation for the effective diffusion is remarkably simple. In Sec. |4] we outline 
three case studies for unstirred mixtures. We demonstrate that modulating the surface can 
increase the reaction yield. We verify that this result is independent of the reaction type 
by obtaining a similar result for the bistable reaction function. Finally, in Sec. |5]we present 
our conclusions. 



2. THEORETICAL FORMULATION 

In this section we describe the mathematical formulation for the autocatalytic reaction. 
We derive the rate equation for homogeneous concentrations, and then generalize this result 
to spatially-varying concentrations on generalized two-dimensional surfaces. We pay close 
attention to understanding the effect on the reaction rate when the surface itself varies with 
time. 

In the homogeneous case, the evolution of n chemical species whose concentrations are 
given by the vector C = (ci(t), C2(t), c„(t)) can be placed in the following canonical form: 

f^^(C). (2.1) 

where the n- vector of functions F is the reaction term describing the interactions between the 
different species. We are interested in particular in the autocatalytic reaction C1 + C2 — > 2c2] 
the dynamical system for such a reaction is given by the equation pair 

where /3 > is the reaction rate. The implied equation d{ci + C2)/dt = is a statement 
of molecular conservation. This system is reduced to a single equation by defining a new 
variable c = C2/ {ci + C2), giving rise to the logistic growth law 

^ = ac(l-c), (2.3) 

where F(c) = ac{l — c) is the reaction function and a = /?(ci + C2) > is the as- 
sociated rate. The evolution of this relative concentration is a contest between linear 
creation and quadratic destruction, which manifests itself through the sigmoid solution 
c = c(0)e'^*/ [1 + c(0)(e'^* — 1)], where c(0) > is the initial concentration. 

The two critical points satisfying dc/dt = are the states c = and c = 1, respectively 
indicating when species C2 is extinct and when the whole space of concentrations is equal 
to that of the product C2; i.e. Ci is extinct. Since dF/dc = a{l — 2c), c = is unstable 
and c = 1 is stable. The phase portrait of the one-dimensional dynamical system is easily 
envisioned: there is a quadratically-increasing reaction away from the repeller c = 0, towards 
larger values of c. When the product is half of the whole concentration mix c = 1/2; the 
reaction rate is maximal, and thereafter it decreasess quadratically towards the attractor at 
c = 1. Hence, if we start with a soup composed of only ci-molecules, the addition of any 
amount of the second species (no matter how small) leads to the annihilation of the first, 
while in the reversed scenario it is what is added, in this case the Ci-molecules, that will be 
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destroyed. Any mixture of the two tends towards a homogeneous soup wholly composed of 
the species C2. 

We can generalize the mass-action law to the inhomogeneous case by using the following 
assumptions: 

• There is diffusion of concentration, arising from thermodynamic fluctuations; 

• There is a large-scale imposed stirring, modelled as an advecting velocity; 

• The substrate Ai on which the substance is placed is modulated as a function of time. 

We write down a continuity equation that takes account of these features. The approach 
we take was discussed by Aris [12j, although the application to reaction-diffusion systems 
is new. We examine the mass balance for the chemical concentration with respect to a 
control area S (t). We work with a general two-dimensional manifold Ai. At time t = the 
manifold is endowed with co-ordinates a^, such that 

M (0) = {xeM:'^\x = x {al, al)}, 

These co-ordinates can be used to label the fluid particles at time t = 0. As time evolves, 
the fluid particles are advected by an imposed flow U, and the particle labels develop time- 
dependence, a{t). We introduce another set of co-ordinates denoted by qit). These are 
fixed in the sense that a point in M^, if located by the co-ordinates q (t), has an instantaneous 
velocity wholly normal to the surface A4 (t). Since the manifold varies smoothly in time, 
there is a set of transformations connecting these co-ordinate systems: 

q = q{a,t), a = a{q,t), ao = a (g (0) , 0) . (2.4) 

A particle that is advected from an initial point ao by the imposed stirring U therefore has 
a velocity 

dq\ _ dq 



Uiq,t) 



where, by Eq. (2.4), the time derivative is taken at fixed particle label a. The last piece of 



formalism needed is the prescription of a metric tensor: 

dx dx 

9ij ^) = ^ ■ ^ = ^(^) ■ 9U)^ 9 = det (gij) , (2.5) 

and the fixedness of the co-ordinate system q is thus made manifest: 

f dx\ dx 



Using the transformations (2.4) and the metric tensor (2.5 ), we obtain a convenient definition 



of area, either as an integral over a fixed domain, or a time-varying one: 

/ dS = \/\g\dq^dq^ = / ^/\g\Jda^da^ , 
Js(t) Js(t) Js(o) 
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where S (0) is the pre-advected domain and 



d ('a^ a"^) 



J 



is the Jacobian of the transformation. This formahsm facihtates the derivation of an analog 
of the Reynolds transport theorem for a concentration field c {q {t) , t) 



d 
dt 



cdS 



S(t) 



S(t) 



Sit) 



dS, 



dc 
dt 



^ + div (C7c) + c log 



dS. 



(2.6) 



This change in the amount of concentration in the control patch must be matched by the 
diffusive flux through the boundary of the patch, and by the amount of matter created or 



AreaS(0 




x(q{r),t) 



FIG. 1: A schematic diagram showing the flux of matter into, and out of a patch of area S. This 
flux comprises advective and diffusive parts, and is given hy J = Uc — k grade, and M is a line 
element on the curve formed by the boundary of the area S. 

destroyed by the reaction (shown schematically in Fig. [T]), that is, 

/ c{q,t)dS = - Kgj:adc-di+ / F (c) dS, 
Js(t) JdS(t) Js(t) 



d 

dt 



where k is the (constant) diffusion coefficient and dS is the boundary of S. A simple 
application of Gauss's law then gives 



d 
dt 



c {q, t) dS 



Sit) 



nAcdS 



Sit) 



F (c) dS, A = div grad. 



(2.7) 



Sit) 



Combining Eqs. (2.6) and (2.7) gives the following local conservation law: 



dc 



^ + div (Uc) = K Ac + F(c) -c- ^ 
at dt 



(2.8) 
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The divergence term div (?7c) can be re- written as Wdq.c = U ■ VgC for incompressible 
flows. We call the term — c(^d log ^/g/dt) the geometric sink: its inclusion is necessary to 
conserve the total number of particles on a time-varying substrate. Note that in previous 
applications [TUl [TT| the scale function was taken to be a growing function of time, and hence 
this extra term was indeed a sink; here we consider a general growth function, and thus this 
term can also act as a source. The geometric source / sink has the following interpretation: 
given a concentration equation for the number of particles per unit volume, a local source 
can be introduced in two ways. The first and more obvious way is to inject particles into 
the system. Here, instead of increasing or decreasing the local number of particles, we 
stretch or squeeze the local area element, so that the local concentration changes. This 
effect vanishes upon integration, so that the total number of particles is conserved in a 
global sense, although locally, the number of particles per unit volume changes because the 
volume itself changes. 



We non-dimensionalize Eq. (2.8) to understand the relative importance of stirring, dif- 



fusion and reaction kinetics on the dynamics. Given a characteristic length scale L, and a 



characteristic speed U of an incompressible velocity field U, Eq. (2.8) is parametrised by 



the Peclet number Pe = LU / k and the Damkohler number Da = aL/U, such that 

^ + div (C7c) = Pe-^Ac + Da c(l - c) - c ^^"^,^ ; (2.9) 
ot at 

these groups respectively describe the ratio of the advective/diffusive and chemical/advective 
timescales. For any stirring variation, the group DaPe = aL'^/n = Const, is unchanged. 
Hece for a fixed diffusion coefficient k, the chemical timescale 1/cr is the only free param- 
eter controlling the dynamics. Thus, by keeping the other parameters fixed, the chemical 
timescale sufficiently represents the dynamical timescale of the reaction-advection-diffusion 
equation [13]. Accordingly, the case of a homogeneous concentration field discussed above 
( 2.3[ ) is equivalent to having a very fast reaction rate with respect to a fixed diffusion of 



order 0{Pe) = 1 and no advection, i.e. Da 3> 1. From now on we shall fix the diffusion 
rate at order unity and describe the dynamics in similar terminology, such that the reaction 
takes place at a rate measured against the diffusion rate. 



3. SCALE SEPARATION WITH STIRRING AND SURFACE MODULATION 

The macroscopic behaviour of a system with phenomena occurring at various length 
and time scales can be described by homogenization theory. The PDE (and its boundary 
conditions) that describes the system, is analyzed as having rapidly oscillating differential 
operators corresponding to the different scales of the phenomena. Taking the appropriate 
limit of infinite scale separation, the solution of the homogenized PDE (known as the cell 
problem) describes the large-scale behavior induced by the small-scale dynamics. In this 
section, we use this approach to calculate the effective diffusivity for the reaction-advection- 
diffusion equation, for a certain class of substrates. The type of substrate modulation we 
specify is rather restrictive, although our results will demonstrate the qualitative effects of 
substrate modulation. 
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3.1. The uniform solution 



If the metric tensor in the reaction-advection-diffusion equation (2.9) satisfies certain 
restrictions, a homogeneous concentration field Co(t) exists. To see this, we set the operators 
div and A to zero in the equation, which gives rise to the following form: 



dco 
~dt 



Da Co (1 - Co) - Co 



dt ■ 



(3.1) 



If, either 



1^1 



p {t) G (g), or more restrictively, 

then we call the metric tensor separable or fully separable respectively. We call the time- 
dependant separable term p{t) the scale function. Note the following observation: 

det{gij) = pG does not imply that gij = pGij. 



One example of such non-equality is furnished by the metric of the torus in Eq. (3.14). By 
choosing an appropriate modulation of the inner and outer toroidal radii, the metric of the 
torus can, however, be made fully separable. Thus, using a separable metric, we have the 
equality 

dlog ^{q,t) ^ 1 dp{t) 

dt p{t) dt ^ ' 

and hence the differential equation for the uniform concentration is itself uniform, and this 
approach is self-consistent. Restricting to this class of substrate modulation, the homoge- 
neous concentration field has an explicit solution as the solution of a Bernoulli equation '■ 

CO {t) = ...^ , • (3-3) 



1 



Da /*{e/o(^-^ 



co(0) ' JO 

In the following applications, we shall make use of the growth function 



7 (t) = Da 



t 



Co (s) ds 







which satisfies the following important result: 



Proposition 1 (The growth function 7 (t) is bounded in time.) This result holds 



when the scale function p (i) is hounded below, < pmin 
perfect derivative, using the notation f{t) := p^^(t)e^"* 



< p{t). We re-write Cq (t) as a 



co(t) 



1 dlog {^^+ Da fl^f{s)ds) 



Da 



dt 



so that its integral is simply 



Co{s)ds 



1 
D^ 

1 
'Da 



log 
log 



+ Da 



Co(0) 
1 + Daco (0) 



f{s)ds 
f{s)ds 



- log 



Co(0) 



8 



Thus, 



7(t) = Dat -2lo. 



1 + DacoiO) / fis)ds 
Jo 

2| log (e^"*/2) - log 1 + Daco (0) ^ f{s)ds 
2 log 



1 + Daco (0) Lfis)ds 



< 21o£ 



^Dat/2 



1 + f * e^'^-ds 

Pmin JO 

and i/izs last quantity has a t-independent upper bound, 70. 



3.2. The homogenized solution 

We homogenize the advection-reaction-diffusion equation on a manifold Ai. For this 
approach to work, we must specialize to a manifolds with certain special properties. Because 
homogenization theory is, in general, applicable only to periodic domains, the surface Ai 
must either be a periodic surface embedded in M.^, or the torus T^, with an appropriate 
modulation of the radii. We make the further restriction that the metric tensor of the 
manifold Ai be fully separable, and that the modulation of the manifold is periodic in time. 
The separability condition means that the co-ordinates q are independent of time. Then, 
the advection-reaction-diffusion equation, written in co-ordinate form q = (g^,^^), is the 
following: 



dc Id 



\G\U'c 



Pe'^ d ( I — dc 



F(c)-c 



dlogp 
~dt~' 



Thus, we have isolated the modulation terms, and the equation can be expressed in a form 
where the differential operators are independent of time: 



I + d,v, (U c) = ^Pe- A„c + f (c) - c^. 

ot p [t) at 



(3.4) 



where Ao0 = p(0)"^ \G\-^'^dq^ {\G\^/^G'Wq^(P) . For ease of notation, we absorb the prefactor 
p (0) into the inverse Peclet number. 

Solving the full logistic model (3.4) is problematic, as it is a non-linear equation. Fortu- 



nately, it is possible to understand the distribution of spatial variations in detail simply by 
studying the linearized form 

c{q,t) = Co{t) + 6ij{q,t), (3.5) 

where 5 -C 1 and the advection only affects the second term, ilj{q,t). Then, the reaction- 
advection-diffusion equation becomes 



|-^Ml-2co(t)] + ^|^ 



-divo (Uij) + ^^AoV^. 
pit) 



(3.6) 



This evolution equation provides a uniform bound on WipW"^, and hence the decomposi- 



tion (3.5) is valid for all times. 
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Proposition 2 (The quantity ||^||2 is uniformly bounded.) To see this, multiply 
Eq. (3.6) by ip^y\g\dq^dq'^ and integrate. Using the incompressibility condition, this gives 



the equation 



where 



d_ 
di 



Pe- 



Then, 



Integrating, 



d 
dt 



Pit) 



2^^ I 



|grad^||2 + Da(l-2co (t)) 



2) 



=d'y{t)/dt 



g\dq^dq^. 



-i{t)\ 



IGMgW < 11^112 W < 11^112 (0) e^W < ll^ll^ (0) e' 



since the growth function 7 (t) is hounded. This gives the required result. 



Before outlining the homogenization method, we re- work Eq. (3.6) such that we are left 
with the simplest possible equation. First, by re-defining time, r = j dtp the equation 
to homogenize has time-dependence only in the velocity and reaction terms: 



dc 



-sir 



^= [-V{q,T)-V, + Pe-'Ao] ^, 



where 



s(r) = Dap(r)(l-2co(r))- 



dlogp 
dr ' 



V' = U' {q,T)p{T). 



The time-dependent source term can be eliminated altogether through the use of the equation 



[-V{q,T)-V, + Pe-'Ao] ^, 



(3.7) 



where 



ijj = {q,T) e 



! s(T)dT 



Note that the transformation r = ds/ p (s) is well-defined because it is invertible, dr/dt = 
p {t)~^ > 0. We now focus on homogenizing Eq. (3.7)). To do this, we must distinguish the 
time- and length scales in the problem. We work with a velocity field V {q, t) that varies 
rapidly in space and time, on scales e and e^, respectively. The parameter e is obtained 
from the ratio e = i/L, where e is the correlation length of the velocity field, and L is the 
domain size. We have chosen the timescale of variation to be e'^ for definiteness, although 
other, slower and faster, scales are possible. We assume that p (t) is periodic in time and 
varies on the scale e"^ while and Gij (q) varies only on the small length scale e, or on the 
box-size scale, but not on both scales. 
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The metric has small-scale variations: We introduce the auxiliary function \l/£ 
\1/ (g/e, r/e^), which satisfies the equation 



dr 



[-6-'Viq,T)-\/, + Pe-'Ao] ^e- 



(3.1 



Next, we introduce auxihary independent variables Q = q/e and a = r/e'^, such that 



A = Aq + 2£:"MivQgrad +£:"^A 



while the time derivative now also has two components: 



_d_ 



+ e~ 



,d_ 
da' 



The PDE now reads, 



where 



Co = -r(Q,a)^ + Pe-iAQ, 

A = (g,a)^ + 2Pe-MivQgrad^, 



Pe-'A^. 



We expand the function in powers of e, as 

g, r, a) = ^o(Q, q, r, a) + e^i(Q, g, r, a) + e^^slQ, 9, r, a) + 



Equating powers of e in the expansion of the equation (3.8), we obtain the following triad 
of problems: 



■9^0 
da 

8^2 <9^o 



da 



+ 



dr 



jCo'^o = 0, 

Co^2 = (£1^1 + A^o) • 



(3.9) 
(3.10) 
(3.11) 



By multiplying Eq. (3.9) by ^/\G\ (Q) \I'o {Q, cr) and integrating, we obtain 



1^ 



0II2 



(Oro)-||^o||2(0) 



|gradQ^'o||2'^o-- 



Thus, if \l/o is (To periodic, then \l/o is independent of Q, and hence, is independent of a. In 
symbols, 

^0 = ^0 (g,r). 
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Using this solution, the second equation (3.10) becomes 



da 



-Co^i = -V'iQ,a) 



This is solved by the ansatz 



where 6"^ {Q,a) solves the cell problem 

de 



da 



da / e'{Q,a) y/\G\dQ^dQ' = 0, 
Juq 



and where 6^* {Q,a) is periodic in Q and a, with periods 1 and cxo. Finally, Eq. (3.11) has a 
solution provided 



JQq 



' Or 



+ £1^1 + £2*0 



\G\dQ'dQ' = 



that is, if 



(r,g) = Pe-^A,v]>o (r, g) + 



1 



The second term can be re-arranged as 



da / Cl<^J^^/\G\dQ'dQ\ 



[ " da [ ^/\G\ (Q) dQ'dQ^ 



2Pe-i / d 



\G\ \dQ 



Qgk 







da I ^/\G\iQ)dQ'dQ^[-V'eq^^^o{q,(T) 



dq^ dq 



where the result on the second line is a consequence of incompressibility and the periodic 
boundary conditions. The homogenization matrix 



— ^ r da j ^\dQHQ^ [-Ve^] 
|"q|o"o Jo JQ.O 



J 



is constant because \/\G\ = \/\G\ (Q). At lowest order, the equation for \l/o (g, t) is 



Or -1 



dq^ dqi 



(3.12) 



and thus the full solution (to Eq. (3.6)) i}) has a slowly-varying spatial component, multiplied 
by a rapidly-varying temporal envelope: 



fs{T)dT 



We note the following result: 
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Proposition 3 (Regularity of the solution \I/o (9)''")) The operator K, is negative in the 
following sense: 

dq^dq'^i)}Cip < 0, 



for any square-integrable function ip, and hence, the solution to the equation (3.12) is 
uniformly hounded in space and time. 

The proof of this statement follows from a straightforward computation based on the J - 
operator: 

dq^dq^i) X' — —^ 



dq" dqi 



--J 





O"0 



dQ'dQ' / dq^dq'^\G\ (Q) W (Q) (Q) i) (q) 



d_d_ 
dq^ dq^ 



Jo.. 



dQ'dQ' / dq^dqW\G\ (Q) (Q) (Q) Wi {q) Wj (q) , 



/ / dQ^dQ^ / dq^dq^^/\G\{Q) {U ■ w) {6 ■ w) 

Jo Jnn J^a 



where Wi = dip/dq^ . We introduce the quantity 6^ = w ■ 6 , which satisfies {d^ — Cq) 6^ 
—U ■ w. Then, 



dq^dq'^ipj'ip 



/•(JO p p 

/ da / dq^dq^ / dQ^ dQ^ ^/\G\ (Q) 

Jo JUg JQq 



d_ 



0w ~ dw^O^^w 



— I da I (ig^(ig^||gradn6'ii)||2 < 0. 



Since \l/o is uniformly bounded in space and time (r), our homogenized solution is bounded, 
and thus consistent with the boundedness result of Prop. [2| provided 

Js{T)dT ^ jm/p{t)]dt 



is bounded in time. This is certainly the case, since the growth function 7 (t) is bounded 
(Prop.g. 

The metric has domain-scale variations: As before, the PDE to homogeninze reads: 



where 



Co(P = -divQ (V (Q, or) 0) + Pe-i Aq0, 
£i0 = -div, {V (Q, a) 0) + 2Pe-MivQgradq0, 
= Pe~^A„<p. 
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The difference is that now \G\ = \G\ (q), where q is the macroscopic scale. This forces the 
velocity field to have the behaviour {Q,q,a) = V'^ {Q,cr) /^/\G\{q), by incompress- 
ibility. The triad of problems (3.9)-(3.11) is unchanged, and the solution to Eq. (3.9) is 
again a function only of the macroscopic scales: 



Using this solution, the second equation (3.10) becomes 
9^1 



da 



-dw,{V{Q,q,a) ^0 (g,r)) 



^0 d 



\G\ dq' 
^0 d 



\G\V' 



V 



This is solved by the ansatz 



where 6'* {Q, q, a) solves the cell problem 

/■CO r 

Cq9' = -V\ / da ¥ (Q, g, a) dQ^dQ^ = 0, 
Jo Jnn 



06' 
9^ 



9'{Q,q,a) d^o 



dq^ 



that is, 



86' 86' 
8^^ 



\G\8Q^ 



Pe-^And' = -V' 



The solution 6' must also satisfy the periodicity condition in Q and a, with periods 1 and ao, 
respectively. Note that the macroscopic variable q appears parametrically in the equation 
for 6'. Finally, Eq. (3.11) has a solution provided 



JnQ L 



8t 



+ £1^1 + £2^0 



dQ^dQ^ = 0, 



that is, if 



^iT,q) = Pe-'A,^oir,q) 



r da [ Ci^!i{q,Qv)dQ^dQ\ 
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The second term can be re-arranged as 



— ^ r da [ dQ'dQ^ 
|"q|o"o Jo Jqo 



d 



d 



9^0 09^ 



^2Pe^ / dQ'dQ^^G^^^^ 



[ " da [ dQ^dQ^ ( 
\^''qWo Jo Jqo 



G\ (q) dq^ ^\ (q) dqi " ^\^\ dq^ dq' 

9^ 9^0 " 



d 



\G\ dq' 
1 



d 



1 

W\ 



da dQ^dQ^ I V' 
I^qIc^o Jo Jn. 



G\{q) dQ' ^\{q) dq- 
dq 



dq^ 

1 d J'^ d 



where the vanishing of the term on the second hne is a consequence of the periodic boundary 
conditions in Q. Thus, we have obtained the homogenized matrix 



da / dQ^dQ' -V'9^ (Q, q, a 



I^qWo Jo 

which determines the equation for \l/o {q, r): 



Or 



IC^o. 



\G\ dq^ xAGldq 



(3.13) 



As before, the full solution (to Eq. ( 3.6[ )) ip has a slowly-varying spatial component, multi- 
plied by a rapidly-varying temporal envelope: 

^(g,a,r) = vl/o (g,r)e^^M'^^ 

where the temporal envelope can depend on a and r. The diffusion operator is negative 
because 

[ ^/\G\dq^dq'^ip}Cip < 0, 
Jng 

for any square-integrable function ip. This can be easily verified by examination of the 
operator 

1 d J'^ d 



\G\dq^ ^^\dq^ 



which satisfies the relation 



\G\dqW4^^ ^ 



\G\ dq' ^\dq 



'dq'i) 



dq^ ./\G\ dq^ 



— dq dq 



1 2 J'^ dilj 



\G\ dq^ dqi 
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Using the trick demonstrated in Prop. |3| the negativity of J", and hence /C is estabhshed, 
and thus \l/o is uniformly bounded in space and time (r and t), consistent with Prop. [2] 

In conclusion, we have obtained effective diffusion operators for the concentration \I' in 
limits when the substrate variation is either large or small in spatial extent (Eqs. (3.12) 
and (3.13)). A result for a metric that varies on both scales is also obtainable from a com- 
bination of these two approaches, provided the metric is separable, \G\ = Gi {q/e)G2 {q)- 
We now turn to the calculation of the diffusion constant for a given manifold. 



3.3. Shear flow on a modulating torus 



We choose an isotropically time- varying torus as our manifold = T^, where both 



the outer and inner radii vary with time according to the constraint 



m 

r{t) 



a > 1. The 



strict inequality is to preserve the topological character of the torus as a ring torus, thus 
preventing any degeneration into a horn- or spindle-like surface. Isotropic motion means 
that neither one of the toroidal angle coordinates is time-dependant. Moreover, either or 
both of the toroidal radii can be regarded as controlling the scale function p(t). We introduce 
orthogonal co-ordinates d — (p, such that x ■ z = r siny9, where x is the position vector and 
z is the constant unit vector in the 2;-direction. Thus, the co-ordinate describes changes 
in angle around the minor circle. The line element is then 



ds' = r^dd'^ + {R + r cos d'/ 



and thus the metric is diagonal: 



9ij 



10 


,0 {R + rcos§f 



(3.14) 



Hence, ^J~g = r{R + r cos'd) = r'^{a + cos-i?). We introduce the scale function: 

Pit) ■■= ritf 



. _ Ritf 



(3.15) 



which is determined by the time variation of either radius. Now the Laplacian on a torus 
takes the form 



A 



1 92 



+ 



1 



92 



sin't? 



d 



r2ai?2 {R + rcos^f dif"^ r {R + r cos'd) d'd' 



which after utilizing the relation (3.15) can be simply written as 

1 92 



P P 



+ 



sin'd d 



(3.16) 



We choose a simple shear flow that varies on the fast scale Q = {ed,e(p). Specifically, 

1 



V 



\G\ 



(o,6(gi)). 
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The cell problem to solve is thus 



d0^ 
d0^ 



+ 



+ 



b de^ 

b 89^ 



da ' ^/\G\dQ^ 
The solution is = (0, 9^ (Q^)), where 



- Pe-'AQ9' 

- Pe-'Ao9' 



0, 



-b iQ') . 



D^^^b{Q'). 



Thus, the matrix J^^ is constant and is equal to 



J 





J^22 



J 
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d9^ 


2 




dh 




dQ^ = Pe 


f 


/o 




/o 


dQ^ 



dQ' 



Note that only the trivial solution is possible if b is replaced with p (a) b {Q"^)'- in this case, 
the time-periodicity forces us to choose the zero-solution, and J = Q. For the non-trivial 
case, the homogenized diffusion equation is 

1 + Pej22 Q2 



sint? d 



a + cos 1? 



(a + cos^?)^ 9v?2 

which can be solved using standard techniques for self-adjoint operators on bounded do- 
mains. 



3.4. The extinction of the catalyst 

The reduction of the advcction term provides a method for understanding the extinction 
problem. We want to know under what circumstances a small initial concentration c (g, 0) = 
5^/^(9,0), 5 <C 1 will go extinct. The linearized homogenization theory just developed is 
appropriate here. The homogenized solution is 

where ^'o satisfies a diffusion equation 

Since /C is negative and the manifold M. is smooth, there is a complete set of eigenfunctions 
{A,^ (q)}^0' "^i^^ corresponding eigenvalues {— A^j^g- Thus, the solution is given by the 
superposition 

i^{q,t) = —7w2^CnA^{q)e J°p('\ 

where the C^'s are constant. For a periodic modulation pmin ^ p(^) < Pmax, the amplitude 
decays to zero if 



min«.A? > DaPe 



P(0) 



The parameter range of extinction is modified in two ways: 
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Da-Pe-^^Al 

P max 





FIG. 2: A schematic diagram showing the possibihty for the extinction of the catalyst, when 
min^A^ > -Da Pe [pmax/p (0)]. The time- varying scale factor Pmax ^ p(0) shifts the parameter 
range of extinction to the right, making extinction less likely. On the other hand, the effective 
diffusion typically shifts the magnitude of the eigenvalues to larger values, thus moving the system 
further into, or closer to the domain of extinction. 



The factor p„ 
Fig. [2} 



> p (0) shifts the extinction threshold to a higher value, as shown in 



The contribution of the effective diffusion in Eqs. (3.12) and (3.13) is negative in 



sign, and thus typically increases the magnitude of the eigenvalues relative to the 
unstirred value. The minimum eigenvalue is then shifted rightward, thus promoting 
extinction. 

Having demonstrated how the full problem can be reduced to a diffusion-type problem in a 
variety of different situations, we turn our attention to the numerical simulation of reaction- 
diffusion equations on various surfaces. 



4. NUMERICAL CASE STUDIES 



In this section we examine three cases in which substrate modulation affects the reaction 
propagation. We focus on a modulating torus in three dimensions, where the determinant 
of the metric tensor is separable, and on a standing-wave substrate in two dimensions, 
in which case the determinant of the metric tensor is not separable. We also examine the 
bistable reaction on the torus, in order to verify that our conclusion, namely that appropriate 
substrate modulation enhances the reaction yield, is independent of the details of the reaction 
kinetics. 

We use both analytical and numerical techniques. Our numerical scheme is a a semi- 
implicit spectral method in two dimensions. In both situations, the PDE to solve can be 
written as 

do . ( (9^c dc dc . \ 

Ac + /i, Ai = At 7r^,7r'^'C, . (4.1) 



ot \oy OX oy 

For the torus, we make the identification x = 6, y = (p, while for the substrate the variables 
X and y have their usual meaning. We have also scaled all lengths in the problem according 
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to an appropriate length scale L, and the time scale is taken to be L'^/D. For the toroidal 
case, L = ro, a characteristic radius, while for the substrate, L is the periodic box size. 
Given a time-periodic modulation, there are three non-dimensional frequencies (timescales) 
in the problem: the diffusive frequency, here normalized to unity, the modulation freqeuncy 
tu, and the reaction frequency a = L'^a/D. We are interested in cases where the effects of the 
chemical reaction and the modulation greatly exceed the diffusive effects, and we therefore 
take u; ~ 5" ^ 1. Following standard practice, we henceforth omit ornamentation over non- 



dimensional quantities. We Fourier-transform the PDE (4.1), which takes the semi- implicit 
form 



dc, 
~dt 



r 

^ (tn) = -k'^Ck (tn+l) + l^k (tn) , /ifc (tn) = / (f Xe''^-'^ 



which can be integrated forward in time using standard techniques (a similar approach 
has been used in solving the Cahn-Hilliard equation; see [2Z]-) This method relieves the 
severe constraint on the timestep arising from a fully explicit treatment of the diffusion, and 
promotes numerical stability. Following standard practice, the discretization in space and 
time is refined until convergence is achieved. 



4.1. A separable modulation on the torus 



We use the toroidal co-ordinate system outhned in Sec. |3 3.3[ with the following radial 
modulation: 

R{t) = ^^V^ < £ < 1, 

^ ^ 1 + £ sm {ojt) 

r{t) = ^, a>l, (4.2) 
a 

with a magnitude e and constant angular frequency u, which gives the pulsation protocol 



r, 



Pit) = .V (4.3) 

[1 + esm {cut}] 

It should be noted that the torus, like the sphere, is special in the sense that the scale 
function, which completely determines the geometric sink, presents itself neatly in the form 
of the radius. Thus, this case can easily be translated to a time-varying sphere. The 



protocol (4.3) modifies the yield of the reaction in the homogeneous case. The yield or 
production is the time-average of the number of molecules created in the reaction, and is 
defined by writing down the differential number of molecules in a a patch of area 

dN (t) = Co (t) dA (t) . 

The integral of this quantity is the instantaneous yield: 

N (t) = An^ap (t) cq (t) , 

while the time-average of this quantity is the mean yield: 

1 '•^ 



{N) = hm - / N{t) dt. 

t— >oo _Z 
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This mean yield is readily computable in the small-e limit. To do this, we use the pulsa- 



tion relation (4.2), together with the homogeneous solution (3.3), to obtain the asymptotic 
relation 



p{t) Co (t) 



1 + [5 sin {ujt) - cos (cut)] + \ [(t2 + Auj^ - cos {2ujt) - 2uja sin {2ujt)] ' 



as t ^ oo. (4.4) 



This gives rise to the long-time average 



(iV) = Afihla hm i / 



N{t)dt, 



An'^ria lim — 



dt 



1 + 



sin (cut) — cos (cut)) + ^^ 



1 - 



cr-^ cos(2a;t)+2a;(T sin(2c<;t) 



47r^roa 



1 ,3^2 -0^2 

1 + -e 



+ 0(e=^), cu>0. (4.5) 



Thus, the amount of product created can either be raised or lowered, depending on the 
reaction rate and the pulsation frequency; for large pulsations, the yield is lowered. The 
exact form of the yield (A^) (cu) is plotted in Fig. |4j Note from this result that (A^) (cu) is 
not continuous at cu = 0. Setting cu to zero and then averaging gives (A^) = 47r^rQa, which 
is the case of no pulsation. On the other hand, for very slow pulsation (compared to the 
reaction rate), Eq. (4.5) becomes 



(A^) = 47r>^a lim - 



T 



dt 



1 + 2e sin (cut) + e"^ sin^ (cut) 



^/l^ + 2e^^ 

1 -£2 



(4.6) 



Equation (4.6) is obtained by setting setting the terms in Eq. (4.5) to zero where cu appears 
as a factor. 

The problem of front propagation on two-dimensional static manifolds has been addressed 
by Gridndrod and Gomatam [28j , and thus some qualitative details of front propagation are 
available. Such analysis involves the reduction of an equation in the Laplace-Beltrami 
operator to an associated equation describing front propagation on the line. Since the line 
is infinite in extent, while the manifold in question is compact, this analysis is valid only 
for intermediate times, after the front has been established, but before the frontal region 
experiences the finite extent of the manifold. We apply this technique to the torus by 
considering an initial disturbance c (■(?, y?, t = 0) = cq {'&) centred on ■(9 = vr. This disturbance 
spreads only in the d direction. Thus c ~ 1 at ^9 ~ vr, while for early times, the reaction 
has not propagated around the torus, and c ~ at ~ 0. We denote the location of the 
front by (t), and introduce a moving coordinate r]{'d,t) = "i? — '^f (t), where the function 

(t) is to be determined. The profile of the front is given by a one- dimensional function: 
c(^?,t) = fiv), and 

dc I d-di 

-f iv) 



dt 



dt 
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Now for < 0, / ^ 1, while for 77 > 0, / ~ 0. Thus, we are interested in a region 77 
where the profile of the function / changes rapidly. We therefore write down the Laplacian 
in the neighbourhood of this point: 



Ac 



sin d dc 



d'd'^ a + cos 7? 

sin di 



R/r 



f" iv) - 



a + cos T^f 

Putting the reaction-diffusion equation together, we have 



/' iv) + o iv) 



f" iv) + f iv) 



ddf sin di 
dt a + cos di 



afil-f) = 0. 



If we stipulate the constant frontal velocity 



k 



d-df 
~dt 



sin 



Const., 



a + cos d{ 

then we are reduced to the reaction-diffusion equation on the line, for the variable rj: 

fr,, + kf, + afil-f) = 0. 



(4.7) 



Eq. (|4.7|) implies that that the velocity of the front is non-constant, and evolves according 

sin-i^f 



to the differential equation 



d^i 
Itt 



a + cos 'di 



{A.i 



In addition to the constant term fc, Eq. (4.8) possesses a curvature-related term that can 



speed up or slow down the front propagation. In particular, there is the possibility of a 
stationary front when 

sin -df 

k + = 0. 

a + cos V{ 

There is no analogue of the steady-state front in react ion- diffusion on the line. Here it 
corresponds to a balance between the tendency of the reaction to propagate, and the curva- 
ture of the torus, which inhibits the reaction propagation. For moderate to large a-values 
a = 10-1000, the curvature term, being a diffusive contribution, is unimportant relative to 
the reaction term, and the dynamical equation for 6f gives approximately linear growth in 
time (this is verified by numerical simulation below). Moreover, in this parameter range, 
it is possible to understand the modulated solution by reference to a fiat-space model. 
Switching on the pulsation clearly will modify the front solution, since the geometric sink 
— c ((9 log y^/c?t) = —cp' (t) /pit) in the concentration equation breaks the Galilean invari- 
ance. However, some quantitative understanding of the front propagation is still possible by 
studying the small-e equation 



dc 
dt 



[l + ed {t)f Ac + (Tc{l 



2ecd' (t) 
[1 + ed (t)] 



2 ' 



(4.9) 



where 

p=[l + ed{t)y^ , d{t) = sin {ujt) , 
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and 

Using regular perturbation theory, it can be shown (Appendix A) that the location of the 
front in this case is given by the formula 

it) = kt + s[A cos {ut) + B sin {ut)] + O (e^) , 

where A and B are amplitudes that are determined from the zeroth and first-order solutions 



of the equation (4.9 ), and kt is the reaction front when e = 0. Note that Mendez ^29j tackles a 
similar problem, but with slowly-varying inhomogeneities and in flat space; the application 
we have in mind has rapidly-varying temporal co-efficients. Thus, the time-evolution for 
a periodic geometric sink is a secular drift, coupled with a local-in-time back-and-forth 
oscillation as the function p {t) is modulated. We can therefore give a qualitative description 
of the front propagation on the modulating torus: there is a secular drift, which is raised or 
lowered over the fiat case due to curvature effects, while there is a back-and-forth oscillation 
in the frontal position due to the modulation of the toroidal area. We turn to numerical 
simulations to check this prediction. 

A numerical approach enables us to describe front propagation in the presence of pulsa- 



tion, and to verify the yield equation (4.5). We work with the pulsation protocol (4.3) and 



choose a ring-shaped disturbance as an initial condition: 

c (t9, ^,t = 0) = e-('^— )'/2'«^ ^ w = 0.2. 

Figure |3] shows the front propagation on the pulsating torus. The catalyst is initially centred 
on "i? = TT and propagates in both directions towards ^9 = 0. The concentration of catalyst 
tends to a uniform amount; however, as the toroidal radii pulsate, the concentration level 
fluctuates. The most striking effect of the pulsating substrate is seen when the yield of 
the reaction is studied, as in Fig. |4} The yield fluctuates over time. The maximum yield 
exceeds the yield in the non-pulsating case, while the minium yield lies below this steady 
value. Fig. |4] (b) shows the time-average yield as a function of pulsation frequency. There 



is a discontinuity at = 0, as discussed in the context of Eqs. (4.5) and (4.6). At slow 
modulation frequencies, the average yield exceeds that of the non-pulsating case, while 
for faster modulation frequencies, the average yield decreases relative to this steady value. 



Fig. |4] (b) also provides a verification of the yield formula Eq. (4.5) and demonstrates the 
concention that surface modulation can enhance the yield. 



In the absence of pulsation, the speed of the front propagation satisfies Eq. (4.8), as 
confirmed in Fig. [s] (a). When the pulsation is switched on, there is still a net drift in the 
location of the front, although locally in time, the front moves forwards and backwards as the 
surface is modulated. This confirmation of our earlier prediction is shown in Fig. [5](b). Next, 
we turn to the study of a qualitatively different case, that of a standing-wave modulation 
on a substrate, in which case the determinant \g\ {x,y,t) is non- separable, and thus the 
geometric sink depends on space and time. 



4.2. A non-seprable modulation: standing wave on a substrate 



In this section we work with a general periodic surface embedded in M.^ with Cartesian 
coordinates {x, y, z). The position vector a; of a point P {x) on the substrate is given in the 
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FIG. 3: Front propagation on the torus at times t = 0, gT, |T, |T, |T, T, where T = 27r/a; 
is the period of the pulsation. The catalyst is initially centred on -i? = vr and propagates in both 
directions towards t? = 0. The concentration of catalyst tends to a uniform amount; however, as 
the torus radii pulsate, the concentration level fluctuates, as shown in subfigures (d)-(f). We have 
taken u! = a = W, and e = 0.5. 



Monge parametrization as 

X = {x,y,f{x,y,t)) , 

where / is a differentiable function of the planar coordinates (x^ = x, = y) and time. The 
metric tensor is thus 

1 + /.' Uy \ 

fjy 1 + /^;' 



^ ~^ fy fxfy 

y 



-Uy 1 + fl 



{9ij) = 

with inverse 

Both of these matrices have determinant 

\9\ = l + f,+fy:=l + {Vjf. 

For a standing-wave surface 

/ (x, y^t) = e sin {kx) sin [ut) , 
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FIG. 4: (a) The instantaneous yield (A^) (t) for uj = 10 = a, and e = 0.5. The system settles down 
to a periodic state wherein the concentration fluctuates homogeneously. The dashed line indicates 
the yield in the absence of pulsation; (b) The time-averaged yield as a function of the pulsation 
frequency. The graph attains its maximum as a; — > 0. However, there is a discontinuity at u; = 0, 
and the yield at zero frequency differs from that for very slow pulsations a; — > 0. This can be seen 



from Eqs. (4.5) and (4.6) At large values of lo, the yield asymptotes to a constant value. There is 



excellent agreement between the yield values provided by this graph and the numerical solution of 
the PDE, and thus, the latter are not shown. 




FIG. 5: Front propagation on the torus, with a = 500. Figure (a) shows front propagation 
on an unmodulated torus, and a comparison with the the front-tracking formula d-d/dt = k + 
[sin'!9/(a -|- cos 1?)], with k = 43. For large times, the comparison is spoiled since the front wraps 
around the torus. Subfigure (b) shows the modulated front. The front drifts with the same velocity 
as in the unmodulated case, while there is a backwards-and-forwards motion as the toroidal area 
varies periodically. 




FIG. 6: Front propagation on the substrate at times t = 0, gT, |r, |r, |r, T, where T = 2t:/u> is 
the period of the pulsation, and uj = a = 10. The catalyst is initially centred on y = and initially 
in the y-direction until a time-periodic state is reached. In reaching this state, the direction of the 
variation changes, as evidenced by subfigures (e) and (f). 



where k and uj are constants, the Laplacian is 



A = — + Cl + — + M^^L 



that is, 



(9^ r , , (9^ e^A;^ sin (/cx) COS (/ca;) sin^ (a;t) d 

^ = TT^ + 1 + ^ ^ cos {kx) sm - ] 2n\ ■ 2 / ^^ V - 

ax"' oy-^ 1 + e^k^ cos^ (kx) sm^ (cjt) ox 

The chemical equation is thus given by 



xt 

dt dx^ ' "-^ ' -""^ dy^ ' 1 + f^dx ' ^ "-^ > 1 + T? 



+ (1 + 7^ + ^^7^ + a (c - c ) - — — c (4.10) 



X 



Figure [6] shows the case of front propagation for an initial concentration level (Gaussian), 
with a wavenumber perpendicular to that of the substrate modulation. The front propa- 
gates into regions of zero concentration, in an inhomogeneous fashion (since there is spatial 
modulation in both directions). After about one period of substrate modulation, the spa- 
tial variation of the concentration field switches from being in the ?/-direction, to being in 
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(a) 



(b) 





(c) 



(d) 



FIG. 7: Front propagation on the substrate at times t = 3T, ^T, and ^T, and ^T, where 
T = juj is the period of the pulsation, w = cr = 10, /c = 27r/L, and e = \. The concentration has 
reached a time-periodic state where ah spatial variations are in the same direction as the direction 
of the substrate modulation. 



the x-direction, aligned with the substrate modulation. Eventually, the system attains a 
time-periodic state, shown in Fig. [7], where the dynamics are driven entirely by the determi- 
nantal function \g\ (x, t\ On the other hand, for front propagation for an initial disturbance 
whose wavenumber is parallel to that of the substrate modulation, the f propagates into 
regions of zero concentration in a homogeneous fashion, and the system rapidly reaches the 
time-periodic state shown in Fig. [7| 

The mean yield is always that associated with with the time-periodic state, since any 
initial configuration tends asymptotically to this state. The yield function is 



We obtain the yield function as a function of the parameters cu, fc, and e by solving Eq. (4.10 ) 



numerically in one dimension [dy = 0). The results are shown in Fig. |8j As before, the mean 
yield as a function of time varies in phase with the substrate modulation, and the maximum 
mean yield exceeds the stationary case. The time-averaged mean yield depends on the 
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FIG. 8: (a) The instantaneous yield {N) (t) for to = 10 = a, k = 27r /uj, and e = 1. This quantity 
is given by the integral J dx J dy^JX + /^c (x, y, t). The system settles down to a periodic state 
wherein the concentration flucuates inhomogeneously. The dashed line indicates the yield in the 
absence of pulsation; (b) The average yield as a function of uj and wavenumber A;, where feo = 2it/L 
is the fundamental wavenumber. In general, for low pulsation frequencies, the yield is raised 
relative to the non-modulated case, while for fast frequencies, the yield is lowered. The higher the 
wavenumber, the stronger the effect. 

frequency of modulation: the slower the modulation, the greater the yield. Increasing the 
wavenumber of the modulation enhances this effect, as seen in Fig. |8] (b). In contrast with 
the toroidal case, the late-time state is not homogeneous, rather it varies periodically in 
space and time, according to the one-dimensional equation 

dc ^ (Pc fj^dc , _ 2N '^ftfxt 
dt dx^ 1 + dx ^ ' l + ff' 

An inhomogeneous final state is undesirable in applications where a pure state involving 
only the product B is required, and thus a pulsation protocol similar to that on the torus is 
preferable over the substrate modulation presented here. 



4.3. The bistable reaction on the torus 

We demonstrate numerically that the reaction yield can be enhanced for other, more 
complicated mass-action laws, such as the bistable reaction. Here, there are two stable 
states c = 0, and c = 1, and an intermediate, unstable state c = ao, where < ao < 1- We 
study this reaction on the pulsating torus; the relevant equation is 



— = Ac+crc(c- 1) (ao - c)-c — , A=- 

at ot p 



1 9^ sin-i? d 

+ 



5^92 (a + cos 9(^2 a + cos ^9^9 

(4-11) 

where p = r [t) ; for our pulsation protocol this is p ^ = [1 + £sin [ujt)] . 

Using a full two-dimensional numerical simulation, we have verified that an arbitrary 
initial state tends either to the state c = 0, or a uniform oscillatory state. The preferred 
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FIG. 9: Characterization of the bistable reaction. Subfigure (a) gives the parameter regimes in 
which either the zero state c = 0, or the oscihatory state, is selected as the asymptotic state. The 
oscillatory state is preferred at small frequencies, and the critical frequency is reduced at large 
ao-values. Subfigure (b) gives the time-averaged yield as a function of uj for ao = 0.4. The time- 
averaged yield exceeds the stationary (u = 0; dotted line) yield for to < iVc, while for uj > uJc, the 
yield is zero. 



state depends on the pulsation parameters and the unstable level ao- To see the relation 
between these parameters, we studied the uniform equation, obtained by setting = = 
in Eq. (4.11). We fixed e = 0.5 and a = 10 and investigated the state selection as a 
function of u and ao. For each value of ao there is a critical frequency such that above 
that frequency, the zero state is preferred, while below that frequency, an oscillatory state 
is selected. This relationship is shown in Fig. |9] (a). For large values of ao, close to ao = 1, 
the critical frequency is shifted downward, indicating that the zero state is preferred for all 
but the slowest of modulation frequencies. We have investigated the time-averaged mean 
yield as a function of u and fixed ao. Fig. |9] (b) shows this relationship for ao = 0.4. 
For u < Uc (ao = 0.4), the time-averaged mean yield exceeds the stationary value (where 
UJ = 0), while for uj > Uc the mean yield is zero. This result demonstrates that while more 
parameter-tuning is required, it is still possible to obtain a yield above the stationary yield 
simply by an appropriate modulation of the substrate. 



5. CONCLUSIONS 

We developed a mass-balance law for fiow-driven chemical reactions on arbitary, time- 
varying surfaces. The derivation is quite general, and takes into account situations where 
the surface co-ordinates are themselves functions of time. Our mass-balance law possesses a 
geometric source / sink, which modifies the reaction. For isotropic surfaces, where the space- 
and time-dependence of the metric tensor are separable, this geometric term is a function of 
time alone, and a homogeneous solution is possible. This solution is explicit for the logistic 
reaction function, and the dependence of the concentration level on the scale function of the 
metric tensor is thus made manifest. 
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In many situations [10], the surface of modulation is isotropic, and this case therefore 
merits close attention. We developed a theory for describing the effects of flow for this 
class of manifold, and for flow fields with small-scale spatial variations. In such a scenario, 
homogenization theory permits one to calculate the distribution of concentration through an 
effective-diffusion equation. Through surface modulation, the effective diffusion coefficient 
depends on space, although this dependence is eliminated for a class of simple shear flows 
on the torus; similar results for other surfaces are easily envisioned. 

Having demonstrated a method for taking account of flow through the use of an effective 
diffusivity, we focused on numerical simulations of react ion- diffusion equations. By numer- 
ically simulating logistic growth and diffusion on the torus, we demonstrated the existence 
of reaction fronts that drift at a constant velocity, but periodically advance and recede, due 
to surface modulation. We also demonstrated that the time-averaged yield of the reaction 
could be increased by surface modulation. A similar result was found for the bistable growth 
law, although careful tuning of the modulation frequency in relation to the bistable param- 
eter is necessary for selection of the required asymptotic state. For non-isotropic surfaces, 
the yield was increased, although a spatially homogeneous state was impossible to attain. 
In summary, our PDE model and its simplifications provide an insight into the simultaneous 
processes of chemical reactions, stirring, and surface modulation, and should prove helpful 
in optimizing the outcome of chemical reactions on variable domains. 
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APPENDIX A 



In this section we calculate the perturbed speed of front propagation for the following 
equation in fiat space: 



dc , , ,,2 ^ X 2ecd' it) / . . x 



where 
and 



p=[l + ed (t)] % d{t)= sin {ut) , 



dx^ 

The amplitude e is assumed to be small, e <^ 1. We define the front as the locus of points 
X{ (t) in the imi-directional solution c{x,t), for which 

c (xf (t) , t) = some constant = Cq. (A-2) 

When e = 0, the front is located at x = xf''^ {t) = kt, where k is the constant velocity, which 
enters into the equation for the front profile: 

~drf dir] ~ <^o) , t] = x - kt, 
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where (—00) = and (00) = 1. We expand the solution to the perturbed problem in 
powers of e: 

c = Co (x, t) +eci {x, t) + O (e^) . 



=<l>o{v) 



The location of the front must change in order for the constraint (A-2) still to be satisfied: 



By Taylor expansion, 



The first-order equation is 

dci 
~dt 



X 



(1) 



Ci {xf\t^ 



dx I 



(A-3) 



Aci + 



dF 

dc 



ci-2[d (t) Aco + d' (t) Co] . 



CO 



We introduce new variables Ci = 0i (77, t). Thus, 



dt 



£^01 - 2 



dit)^ + d' it) iv) 



where 



£^01 - d (t) bi [r]) - d' (t) 62 iv) 

C - — V— — 

dri^ dri dc ^^^(^)' 



Since d is periodic, we can write d = 3? {6oe *'^*) without loss of generality, and thus 0i {i], t) 



-iwt , 



(77), where 0^; [j]) satisfies the equation 



Thus, 



and hence Eq. (A-3) becomes 



(1)__^^0^^_(OK^ 



In other words, 

Xf (t) = (t) + £ [A cos (cut) + 5 sin (cut)] + O (e^) 
where A and -B are constants, as claimed in Sec. |4j 
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